## Callin Switzer
## 29 Nov 2016
## 6 Dec 2016 Update -- re-digitized pollen and anthers very carefully
## may not need any smoothing
## 8 Dec Update: used cross validation for decide smoothing parameter for smoothing spline
## 8 Feb 2017 Update: Cleaned up code for submission as supplemental info
## 22 June 2017 Update: Made example plots for position/speed/accel



## Kalmia pollen and anther kinematics
### 1. Read in digitized files
### 2. Smooth digitized points, using smoothing spline with tuning parameter from cross validation
### 3. Use smoothed points to calculate velocity and acceleration (normal and tangential)

Install packages

ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c("ggplot2", "scales", "multcomp", "plyr", "car", "lme4", "signal")
ipak(packages)
##  ggplot2   scales multcomp     plyr      car     lme4   signal 
##     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE
# print session info and time the code was run last
Sys.time()
## [1] "2017-06-22 20:11:20 EDT"
sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] signal_0.7-6    lme4_1.1-12     Matrix_1.2-8    car_2.1-4      
##  [5] plyr_1.8.4      multcomp_1.4-6  TH.data_1.0-8   MASS_7.3-45    
##  [9] survival_2.40-1 mvtnorm_1.0-5   scales_0.4.1    ggplot2_2.2.1  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.9        nloptr_1.0.4       tools_3.3.2       
##  [4] digest_0.6.12      evaluate_0.10      tibble_1.2        
##  [7] gtable_0.2.0       nlme_3.1-130       lattice_0.20-34   
## [10] mgcv_1.8-16        yaml_2.1.14        parallel_3.3.2    
## [13] SparseM_1.74       stringr_1.1.0      knitr_1.15.1      
## [16] MatrixModels_0.4-1 rprojroot_1.2      grid_3.3.2        
## [19] nnet_7.3-12        rmarkdown_1.3      minqa_1.2.4       
## [22] magrittr_1.5       backports_1.0.5    codetools_0.2-15  
## [25] htmltools_0.3.5    splines_3.3.2      assertthat_0.1    
## [28] pbkrtest_0.4-6     colorspace_1.3-2   quantreg_5.29     
## [31] sandwich_2.3-4     stringi_1.1.2      lazyeval_0.2.0    
## [34] munsell_0.4.3      zoo_1.7-14

Read in data

# set directory that contains datasets
dataDirectory = '/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/DatasetsSupplemental/'

# read in metadata
dfile <- "LaurelsOnly.csv"
metDat <- read.csv(paste0(dataDirectory, dfile))

metDat <- metDat[metDat$redigitizedFile!= "", ]

# set constants:
fps <- 5000 # frames per second

Read in digitized files (xypts)

Smooth xypts

Calculate pollen and anther acceleration and velocity

Note: The smoothing parameter was chosen with 10-fold CV

# read in each .csv file for analysis
digdirect <- "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/DatasetsSupplemental/CleanKalmiaDigitized/"

newDF <- data.frame()



for(ii in 1:nrow(metDat)){
     
     ddfile <- paste0(digdirect, metDat$redigitizedFile[ii])
     dp <- read.csv(ddfile)
     
     # calibrate locations, based on digitized pin or other object
     # calibration points
     pin <- data.frame(dp$pt1_cam1_X, dp$pt1_cam1_Y, dp$pt2_cam1_X, dp$pt2_cam1_Y)
     pin <- pin[complete.cases(pin), ]
     
     # get the number of pixels in the calibration
     PixInPin <- (sqrt((pin$dp.pt1_cam1_X - pin$dp.pt2_cam1_X)^2 +
                            (pin$dp.pt1_cam1_Y-pin$dp.pt2_cam1_Y)^2)) / 
          metDat$CalSizeMM[ii] # to get to mm
     
     
     # get anther and pollen locations
     antherPoll <- data.frame(anthx = dp$pt3_cam1_X, anthy= dp$pt3_cam1_Y, 
                              polx = dp$pt4_cam1_X, poly= dp$pt4_cam1_Y)
     
     
     # get frame where pollen starts and leaves
     antherPoll$polStart = 1:nrow(antherPoll) == metDat$framePollenStartedLeaving[ii]
     antherPoll$polEnd = 1:nrow(antherPoll) == metDat$framePollenReleaseComplete[ii]
     
     # gives only rows where either anth and pollen are complete
     antherPoll = antherPoll[ complete.cases(antherPoll[c('anthx')]) | 
                                   complete.cases(antherPoll[c('polx')]), ]
     
     # if x value starts to right of screen, reverse points,
     # so all x values start on the left part of the screen at 0
     if(lm(antherPoll[,1] ~ I(1:length( antherPoll[,1])))$coefficients[2] < 0 ){
          antherPoll$anthx <- metDat[ii,'vidWidth'] - antherPoll$anthx
          antherPoll$polx <- metDat[ii,'vidWidth'] - antherPoll$polx
     }
     
     
     # cbind data frame, to add smoothed columns
     antherPoll <- data.frame(cbind(antherPoll, antherPoll))
     
     # smooth with a smoothing spline, spar = 0.29 -- this smoothing
     # spline was chosen through cross validation
     foo = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(y){
          yy = na.omit(antherPoll[, y])
          xx = 1:length(yy)
          
          sm1 <- smooth.spline(x = xx, y = yy, spar = 0.29)
          antherPoll[, y][complete.cases(antherPoll[, y])] <<- sm1$y
     })
     
     # add time to data frame
     antherPoll$tme = 0: (nrow(antherPoll) - 1) / fps # time
     
     # add columns with absolute position into dataframe (calculated from smoothed data)
     # calculate position from starting point, not from minimum point
     bar = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(x){
          newName = paste0(x, ".abs")
          tmp <- antherPoll[,x] / PixInPin / 1000
          antherPoll[,newName] <<- tmp - na.omit(tmp)[1]
          #antherPoll[,newName] <<- tmp - min(na.omit(tmp))
     })
     
     # add columns to show velocity, based on smoothed, absolute position
     # velocity is in m/s
     bat = sapply(X = c("anthx.1.abs", "anthy.1.abs", "polx.1.abs", "poly.1.abs"), 
                  FUN = function(x){
                       newName = paste0(x, ".vel")
                       tmp <-  c(NaN, diff(antherPoll[,x])) * fps # add a NaN to beginning of data
                       antherPoll[,newName] <<- tmp 
                  })
     
     # calculate speed
     antherPoll$anthspeed = sqrt(antherPoll$anthx.1.abs.vel^2 + antherPoll$anthy.1.abs.vel^2)
     antherPoll$polspeed = sqrt(antherPoll$polx.1.abs.vel^2 + antherPoll$poly.1.abs.vel^2)
     
     ###########################################
     # pollen acceleration
     polVelocity = cbind(antherPoll$polx.1.abs.vel, antherPoll$poly.1.abs.vel)
     polSpeed = antherPoll$polspeed
     # plot(polSpeed)
     tme = antherPoll$tme
     polAccel = data.frame(rbind(c(NA, NA), apply(polVelocity, MARGIN = 2, FUN = diff))) * fps
     
     # calculate unit tangent vector
     T_t = polVelocity / polSpeed
     
     DT = data.frame(rbind(c(NA, NA), apply(T_t, MARGIN = 2, FUN = diff))) * fps
     NormDT = sqrt(DT[,1]^2 + DT[,2]^2)
     Curvature = NormDT / polSpeed
     
     # compute a_N (normal acceleration) and a_T (tangential acceleration)
     # a_T = ds/dt
     a_T =  c(NA, diff(polSpeed) * fps)
     N_t = data.frame(t(sapply(1:nrow(DT), FUN = function(x) unlist(DT[x, ] / NormDT[x]))))
     # plot(a_T, type = "l", ylim = c(-3000, 3000))
     
     # a_N = speed^2 * curvature
     a_N = polSpeed^2 * Curvature
     
     
     # check total accel by adding normal and tangential accelerations
     # a_total = a_T * T_t + a_N * N_t
     a_total = as.data.frame(t(sapply(X = 1:nrow(polAccel), FUN  = function(x) a_T[x] * T_t[x, ] + a_N[x] * N_t[x,] )))

     a_T_Pol = a_T
     
     # calculate magnitude of acceleration, using two methods (should be the same)
     # 1. Normal and tangential acceleration
     a_mag1 = sqrt(a_T^2 + a_N^2)
     amag2 = sqrt(polAccel[,1]^2 + polAccel[,2]^2)

     
     ###########################################
     # anther acceleration
     anthVelocity = cbind(antherPoll$anthx.1.abs.vel, antherPoll$anthy.1.abs.vel)
     anthSpeed = antherPoll$anthspeed
     tme = antherPoll$tme
     anthAccel = data.frame(rbind(c(NA, NA), apply(anthVelocity, MARGIN = 2, FUN = diff))) * fps
     
     # unit tangent vector
     T_t = anthVelocity / anthSpeed
     
     DT = data.frame(rbind(c(NA, NA), apply(T_t, MARGIN = 2, FUN = diff))) * fps
     NormDT = sqrt(DT[,1]^2 + DT[,2]^2)
     Curvature = NormDT / anthSpeed
     
     # compute a_N (normal acceleration) and a_T (tangential acceleration)
     # a_T = ds/dt
     a_T =  c(NA, diff(anthSpeed) * fps)
     a_T_anth = a_T
     N_t = data.frame(t(sapply(1:nrow(DT), FUN = function(x) unlist(DT[x, ] / NormDT[x]))))
     # plot(a_T, type = "l", ylim = c(-3000, 3000))
     
     # a_N = speed^2 * curvature
     a_N = anthSpeed^2 * Curvature
     
     
     # check total accel by adding normal and tangential accelerations
     # a_total = a_T * T_t + a_N * N_t
     a_total = as.data.frame(t(sapply(X = 1:nrow(anthAccel), FUN  = function(x) a_T[x] * T_t[x, ] + a_N[x] * N_t[x,] )))
     
     # align frames
     tmeRoll <- seq(from = -which(antherPoll$polStart) + 1, length.out = length(tme)) / fps
     
     # update for different visualization
     # align frames -- center, based on max speed
     # tmeRoll <- seq(from = -which.max(antherPoll$polspeed) + 1, length.out = length(tme)) / fps
     
     dfi <- data.frame(anthSpeed, polSpeed, a_T_anth, a_T_Pol, tme, 
                       trial = metDat$VideoName[ii], 
                       tmeStart = antherPoll$polStart,
                       tmeEnd = antherPoll$polEnd, 
                       centeredTime = tmeRoll)
     
     newDF <- rbind(newDF,dfi) 
     print(ii)
     
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32

Visualize pollen and anther speed

theme_set(theme_classic())

savePath = "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/FigsNotForSubmission/"

# anther speed
ggplot(newDF, aes(x = centeredTime, y = anthSpeed, group = trial)) + 
     geom_line(alpha = 0.5) +
     xlim(c(-0.01, 0.02)) + 
     ylim(c(0,6)) + 
     labs(x = "Time (s)", y = "Anther speed (m/s)") + 
ggsave(paste0(savePath, "antherSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).

## Warning: Removed 113 rows containing missing values (geom_path).

# pollen speed
ggplot(newDF, aes(x = centeredTime, y = polSpeed, group = trial)) + 
     geom_line(alpha = 0.5) + 
     xlim(c(-0.01, 0.02)) +
     ylim(c(0,6)) +  
     labs(x = "Time (s)", y = "Pollen speed (m/s)")
## Warning: Removed 53 rows containing missing values (geom_path).

ggsave(paste0(savePath, "pollenSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 53 rows containing missing values (geom_path).
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_anth, group = trial)) + 
     geom_line(alpha = 0.5) + 
     #ylim(c(-2500, 4000)) +
     xlim(c(-0.01, 0.02)) + 
     labs(x = "Time (s)", y = "Anther tangential acceleration  (m/s/s)")
## Warning: Removed 145 rows containing missing values (geom_path).

ggsave(paste0(savePath, "antherTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 145 rows containing missing values (geom_path).
# pollen tangential acceleration
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_Pol, group = trial)) + 
     geom_line(alpha = 0.5) + 
     #ylim(c(-2500, 4000)) +
     xlim(c(-0.01, 0.02)) + 
     labs(x = "Time (s)", y = "Pollen tangential acceleration  (m/s/s)")
## Warning: Removed 85 rows containing missing values (geom_path).

ggsave(paste0(savePath, "PollenTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 85 rows containing missing values (geom_path).

Find max for each measurement for each trial and visualize

# anther speed
mmx <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
     tmp <- newDF[newDF$trial == x, ]
     return (unlist(tmp[which.max(tmp$anthSpeed),]))
})))
mmx$trial <- row.names(mmx)


ggplot() + 
     geom_line(data = newDF, aes(x = centeredTime, y = anthSpeed, group = as.factor(trial)), alpha = 0.5) +
     xlim(c(-0.01, 0.02)) + 
     ylim(c(0,6)) + 
     labs(x = "Time (s)", y = "Anther speed (m/s)") + 
     geom_point(data = mmx, aes(x = centeredTime, y = anthSpeed), color = 'red', alpha = 0.5) + 
     theme(legend.position = "none") 
## Warning: Removed 113 rows containing missing values (geom_path).

#+  facet_wrap(~ trial)
ggsave(paste0(savePath, "antherSpeedMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).
# pollen speed
mmp <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
     tmp <- newDF[newDF$trial == x, ]
     tmp <- tmp[abs(tmp$centeredTime) < 0.01, ]
     return (unlist(tmp[which.max(tmp$polSpeed),]))
})))

mmp$trial <- row.names(mmp)

# pollen speed
ggplot() + 
     geom_line(data = newDF, aes(x = centeredTime, y = polSpeed, group = trial), alpha = 0.5) + 
     xlim(c(-0.01, 0.02)) + 
     ylim(c(0,6)) + 
     labs(x = "Time (s)", y = "Pollen speed (m/s)") + 
     geom_point(data = mmp, aes(x = centeredTime, y = polSpeed), color = 'red', alpha = 0.5) + 
     theme(legend.position = "none") 
## Warning: Removed 53 rows containing missing values (geom_path).

ggsave(paste0(savePath, "pollenSpeedMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 53 rows containing missing values (geom_path).
# anther acceleration
mma <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
     tmp <- newDF[newDF$trial == x, ]
     
     # get only points that are within 0.05 seconds of the centered time
     # to ignore the anthers hitting the other side of the flower
     #tmp <- tmp[abs(tmp$centeredTime) < 0.002, ]
     return (unlist(tmp[which.max(tmp$a_T_anth),]))
})))

mma$trial <- row.names(mma)

ggplot() + 
     geom_line(data = newDF, aes(x = centeredTime, y = a_T_anth, group = trial), alpha = 0.5) + 
     #ylim(c(-2500, 4000)) +
     xlim(c(-0.01, 0.02)) + 
     labs(x = "Time (s)", y = "Anther tangential acceleration  (m/s/s)") + 
     geom_point(data = mma, aes(x = centeredTime, y = a_T_anth), color = 'red', alpha = 0.5)
## Warning: Removed 145 rows containing missing values (geom_path).

ggsave(paste0(savePath, "antherTangAccelMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 145 rows containing missing values (geom_path).
# pollen acceleration
mmpp <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
     tmp <- newDF[newDF$trial == x, ]
     #tmp <- tmp[abs(tmp$centeredTime) < 0.003, ]
     return (unlist(tmp[which.max(tmp$a_T_Pol),]))
})))

mmpp$trial <- row.names(mmpp)

ggplot() + 
     geom_line(data = newDF, aes(x = centeredTime, y = a_T_Pol, group = trial), alpha = 0.5) + 
     #ylim(c(-2500, 4000)) +
     xlim(c(-0.01, 0.02)) + 
     labs(x = "Time (s)", y = "Pollen tangential acceleration  (m/s/s)") + 
     geom_point(data = mmpp, aes(x = centeredTime, y = a_T_Pol), color = 'red', alpha = 0.5)
## Warning: Removed 85 rows containing missing values (geom_path).

ggsave(paste0(savePath, "pollenTangAccelMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 85 rows containing missing values (geom_path).

Make example figure to show position/speed/acceleration

ii = 28
print(ii)
## [1] 28
ddfile <- paste0(digdirect, metDat$redigitizedFile[ii])

ddfile <- "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/DatasetsSupplemental/CleanKalmiaDigitized/20150608_135726xypts.csv"
dp <- read.csv(ddfile)

# calibrate locations, based on digitized pin or other object
# calibration points
pin <- data.frame(dp$pt1_cam1_X, dp$pt1_cam1_Y, dp$pt2_cam1_X, dp$pt2_cam1_Y)
pin <- pin[complete.cases(pin), ]

# get the number of pixels in the calibration
PixInPin <- (sqrt((pin$dp.pt1_cam1_X - pin$dp.pt2_cam1_X)^2 +
                       (pin$dp.pt1_cam1_Y-pin$dp.pt2_cam1_Y)^2)) / 
     metDat$CalSizeMM[ii] # to get to mm


# get anther and pollen locations
antherPoll <- data.frame(anthx = dp$pt3_cam1_X, anthy= dp$pt3_cam1_Y, 
                         polx = dp$pt4_cam1_X, poly= dp$pt4_cam1_Y)




# get frame where pollen starts and leaves
antherPoll$polStart = 1:nrow(antherPoll) == metDat$framePollenStartedLeaving[ii]
antherPoll$polEnd = 1:nrow(antherPoll) == metDat$framePollenReleaseComplete[ii]

# gives only rows where either anth and pollen are complete
antherPoll = antherPoll[ complete.cases(antherPoll[c('anthx')]) | 
                              complete.cases(antherPoll[c('polx')]), ]



# if x value starts to right of screen, reverse points,
# so all x values start on the left part of the screen at 0
if(lm(antherPoll[,1] ~ I(1:length( antherPoll[,1])))$coefficients[2] < 0 ){
     antherPoll$anthx <- metDat[ii,'vidWidth'] - antherPoll$anthx
     antherPoll$polx <- metDat[ii,'vidWidth'] - antherPoll$polx
}



plot(antherPoll$polx, antherPoll$poly, asp = 1)

# cbind data frame, to add smoothed columns
antherPoll <- data.frame(cbind(antherPoll, antherPoll))

# smooth with a smoothing spline, spar = 0.29 -- this smoothing
# spline was chosen through cross validation
foo = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(y){
     yy = na.omit(antherPoll[, y])
     xx = 1:length(yy)
     
     sm1 <- smooth.spline(x = xx, y = yy, spar = 0.29)
     antherPoll[, y][complete.cases(antherPoll[, y])] <<- sm1$y
})

plot(antherPoll$anthx.1, antherPoll$anthy.1)

antherPoll <- antherPoll[10:nrow(antherPoll), ]

# add time to data frame
antherPoll$tme = 0: (nrow(antherPoll) - 1) / fps # time

# add columns with absolute position into dataframe (calculated from smoothed data)
# calculate position from starting point, not from minimum point
bar = sapply(X = c("anthx.1", "anthy.1", "polx.1", "poly.1"), FUN = function(x){
     newName = paste0(x, ".abs")
     tmp <- antherPoll[,x] / PixInPin / 1000
     antherPoll[,newName] <<- tmp - na.omit(tmp)[1]
     #antherPoll[,newName] <<- tmp - min(na.omit(tmp))
})

plot(antherPoll$anthx.1.abs, antherPoll$anthy.1.abs)

# add columns to show velocity, based on smoothed, absolute position
# velocity is in m/s
bat = sapply(X = c("anthx.1.abs", "anthy.1.abs", "polx.1.abs", "poly.1.abs"), 
             FUN = function(x){
                  newName = paste0(x, ".vel")
                  tmp <-  c(NaN, diff(antherPoll[,x])) * fps # add a NaN to beginning of data
                  antherPoll[,newName] <<- tmp 
             })

# calculate speed
antherPoll$anthspeed = sqrt(antherPoll$anthx.1.abs.vel^2 + antherPoll$anthy.1.abs.vel^2)
antherPoll$polspeed = sqrt(antherPoll$polx.1.abs.vel^2 + antherPoll$poly.1.abs.vel^2)

###########################################
# pollen acceleration
polVelocity = cbind(antherPoll$polx.1.abs.vel, antherPoll$poly.1.abs.vel)
polSpeed = antherPoll$polspeed
# plot(polSpeed)
tme = antherPoll$tme
polAccel = data.frame(rbind(c(NA, NA), apply(polVelocity, MARGIN = 2, FUN = diff))) * fps

# calculate unit tangent vector
T_t = polVelocity / polSpeed

DT = data.frame(rbind(c(NA, NA), apply(T_t, MARGIN = 2, FUN = diff))) * fps
NormDT = sqrt(DT[,1]^2 + DT[,2]^2)
Curvature = NormDT / polSpeed

# compute a_N (normal acceleration) and a_T (tangential acceleration)
# a_T = ds/dt
a_T =  c(NA, diff(polSpeed) * fps)
N_t = data.frame(t(sapply(1:nrow(DT), FUN = function(x) unlist(DT[x, ] / NormDT[x]))))
# plot(a_T, type = "l", ylim = c(-3000, 3000))

# a_N = speed^2 * curvature
a_N = polSpeed^2 * Curvature


# check total accel by adding normal and tangential accelerations
# a_total = a_T * T_t + a_N * N_t
a_total = as.data.frame(t(sapply(X = 1:nrow(polAccel), FUN  = function(x) a_T[x] * T_t[x, ] + a_N[x] * N_t[x,] )))

a_T_Pol = a_T




ggplot(antherPoll, aes(y = poly.1.abs, x = polx.1.abs)) + 
     geom_path(alpha = 0.5)  + 
     coord_fixed() + 
     ylim(0, 0.11)+
     labs(y = "Vertical Position (m)", x = "Horizontal position (m)")

x1 <- ggplot(antherPoll, aes(x = tme, y = polx.1.abs)) + 
     geom_line() + 
     labs(x = "", y = "Pollen\nhorizontal position (m)") + 
     ylim(-0.01, 0.11)+ 
     theme(axis.title.y=element_text(margin=margin(0,20,0,0)))+ 
     theme(plot.margin = unit(c(1,1,-1,0), "cm")) + 
ggsave(filename = paste0(savePath,"/horizontalpos.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)

x2 <- ggplot(antherPoll, aes(x = tme, y = poly.1.abs)) + 
     geom_line() + 
     ylim(-0.01, 0.11)+
     labs(x = "", y = "Pollen\nvertical Position (m)")+ 
     theme(axis.title.y=element_text(margin=margin(0,20,0,0)))+ 
     theme(plot.margin = unit(c(1,1,-1,0), "cm")) + 
ggsave(filename = paste0(savePath,"/verticalpos.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)

x3 <- ggplot(antherPoll, aes(x = tme, y = polspeed)) + 
     geom_line() + 
     labs(x = "", y = "Pollen\nspeed (m/s)")+ 
     theme(plot.margin = unit(c(1,1,-5,0), "cm")) + 
     theme(axis.title.y=element_text(margin=margin(0,20,0,0)))
x3
## Warning: Removed 1 rows containing missing values (geom_path).

ggsave(filename = paste0(savePath,"/speed.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)
## Warning: Removed 1 rows containing missing values (geom_path).
x4 <- ggplot(antherPoll, aes(x = tme, y = a_T_Pol)) + 
     geom_line() + 
     labs(x = "Time (s)", y = "Pollen\ntangential accel. (m/s/s)") + 
     theme(axis.title.y=element_text(margin=margin(0,20,0,0)))+ 
     theme(plot.margin = unit(c(0,0,0,0), "cm")) + 
ggsave(filename = paste0(savePath,"/pollenAccel.png"), width = 5*1.5, height = 1.5*1.5, dpi =200)
## Warning: Removed 2 rows containing missing values (geom_path).
#install.packages("cowplot")
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
plot_grid(x1, x2, x3, x4, align='vh', ncol = 1)
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_path).

ggsave(paste0(savePath,"fourplots.png"), height = 10, width = 10, dpi = 120)

Visualize pollen and anther speed

theme_set(theme_classic())

savePath = "/Users/callinswitzer/Dropbox/ExperSummer2016/Kalmia/Manuscript/FigsNotForSubmission/"

# anther speed
ggplot(newDF, aes(x = centeredTime, y = anthSpeed, group = trial)) + 
     geom_line(alpha = 0.5) +
     xlim(c(-0.01, 0.02)) + 
     ylim(c(0,6)) + 
     labs(x = "Time (s)", y = "Anther speed (m/s)") + 
ggsave(paste0(savePath, "antherSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).

## Warning: Removed 113 rows containing missing values (geom_path).

# pollen speed
ggplot(newDF, aes(x = centeredTime, y = polSpeed, group = trial)) + 
     geom_line(alpha = 0.5) + 
     xlim(c(-0.01, 0.02)) +
     ylim(c(0,6)) +  
     labs(x = "Time (s)", y = "Pollen speed (m/s)")
## Warning: Removed 53 rows containing missing values (geom_path).

ggsave(paste0(savePath, "pollenSpeed01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 53 rows containing missing values (geom_path).
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_anth, group = trial)) + 
     geom_line(alpha = 0.5) + 
     #ylim(c(-2500, 4000)) +
     xlim(c(-0.01, 0.02)) + 
     labs(x = "Time (s)", y = "Anther tangential acceleration  (m/s/s)")
## Warning: Removed 145 rows containing missing values (geom_path).

ggsave(paste0(savePath, "antherTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 145 rows containing missing values (geom_path).
# pollen tangential acceleration
# anther tangential acceleration
ggplot(newDF, aes(x = centeredTime, y = a_T_Pol, group = trial)) + 
     geom_line(alpha = 0.5) + 
     #ylim(c(-2500, 4000)) +
     xlim(c(-0.01, 0.02)) + 
     labs(x = "Time (s)", y = "Pollen tangential acceleration  (m/s/s)")
## Warning: Removed 85 rows containing missing values (geom_path).

ggsave(paste0(savePath, "PollenTangAccel01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 85 rows containing missing values (geom_path).

Find max for each measurement for each trial and visualize

# anther speed
mmx <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
     tmp <- newDF[newDF$trial == x, ]
     return (unlist(tmp[which.max(tmp$anthSpeed),]))
})))
mmx$trial <- row.names(mmx)


ggplot() + 
     geom_line(data = newDF, aes(x = centeredTime, y = anthSpeed, group = as.factor(trial)), alpha = 0.5) +
     xlim(c(-0.01, 0.02)) + 
     ylim(c(0,6)) + 
     labs(x = "Time (s)", y = "Anther speed (m/s)") + 
     geom_point(data = mmx, aes(x = centeredTime, y = anthSpeed), color = 'red', alpha = 0.5) + 
     theme(legend.position = "none") 
## Warning: Removed 113 rows containing missing values (geom_path).

#+  facet_wrap(~ trial)
ggsave(paste0(savePath, "antherSpeedMax01_CVSmoothSpline.pdf"), width = 5, height = 4)
## Warning: Removed 113 rows containing missing values (geom_path).
# pollen speed
mmp <- as.data.frame(t(sapply(unique(as.character(newDF$trial)), FUN = function(x){
     tmp <- newDF[newDF$trial == x, ]
     tmp <- tmp[abs(tmp$centeredTime) < 0.01, ]
     return (unlist(tmp[which.max(tmp$polSpeed),]))
})))

mmp$trial <- row.names(mmp)

Estimate means and CI’s for acceleration, and speed

Conduct model diagnostics to evaluate model quality

md  = merge(x = mmx[, c('trial', 'anthSpeed')], metDat, by.x = "trial", by.y = 
                 "VideoName")
md = merge(x = mmp[, c('trial', 'polSpeed')], md, by = "trial")

md = merge(x = mmpp[, c('trial', 'a_T_Pol')], md, by = "trial")
md = merge(x = mma[, c('trial', 'a_T_anth')], md, by = "trial")


#LMER
hist(md$anthSpeed) #dist is fine

modVelMaxAnth <- lmer(formula = anthSpeed ~  (1|plant/FlowerNumber), data = md)
summary(modVelMaxAnth)
## Linear mixed model fit by REML ['lmerMod']
## Formula: anthSpeed ~ (1 | plant/FlowerNumber)
##    Data: md
## 
## REML criterion at convergence: 63.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.52444 -0.65641  0.00419  0.52344  1.89172 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  FlowerNumber:plant (Intercept) 0.2844   0.5333  
##  plant              (Intercept) 0.4103   0.6406  
##  Residual                       0.1608   0.4010  
## Number of obs: 32, groups:  FlowerNumber:plant, 15; plant, 7
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    3.772      0.296   12.74
confint(modVelMaxAnth)
## Computing profile confidence intervals ...
##                 2.5 %    97.5 %
## .sig01      0.2798029 0.9608723
## .sig02      0.0000000 1.2573804
## .sigma      0.2973171 0.5798682
## (Intercept) 3.1531345 4.3868581
# diagnostics
plot(modVelMaxAnth)

# QQPlot for group-level effects
qqnorm(ranef(modVelMaxAnth)$plant[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modVelMaxAnth)$plant[[1]])

# QQPlot for group-level effects
qqnorm(ranef(modVelMaxAnth)$FlowerNumber[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modVelMaxAnth)$FlowerNumber[[1]]) # 

modVelMaxPol <- lmer(formula = polSpeed ~ (1|plant/FlowerNumber), data = md)
summary(modVelMaxPol) # final model for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: polSpeed ~ (1 | plant/FlowerNumber)
##    Data: md
## 
## REML criterion at convergence: 71.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.03862 -0.49111 -0.07546  0.27208  1.95680 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  FlowerNumber:plant (Intercept) 0.2414   0.4914  
##  plant              (Intercept) 0.1210   0.3479  
##  Residual                       0.3089   0.5558  
## Number of obs: 32, groups:  FlowerNumber:plant, 15; plant, 7
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   3.5019     0.2141   16.36
confint(modVelMaxPol)
## Computing profile confidence intervals ...
##                  2.5 %    97.5 %
## .sig01      0.06228879 0.9128825
## .sig02      0.00000000 0.8681086
## .sigma      0.41375779 0.7940168
## (Intercept) 3.05487009 3.9520846
plot(modVelMaxPol)

# log transformed, b/c distribution is skewed.
hist(md$a_T_Pol)

hist(log(md$a_T_Pol)) #better

modAccMaxPol <- lmer(formula = log(a_T_Pol) ~  (1|plant/FlowerNumber), data = md)
summary(modAccMaxPol) #final model for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(a_T_Pol) ~ (1 | plant/FlowerNumber)
##    Data: md
## 
## REML criterion at convergence: 13.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.62436 -0.66206  0.03069  0.61934  2.08762 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  FlowerNumber:plant (Intercept) 0.001824 0.0427  
##  plant              (Intercept) 0.066778 0.2584  
##  Residual                       0.056679 0.2381  
## Number of obs: 32, groups:  FlowerNumber:plant, 15; plant, 7
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   8.3302     0.1095    76.1
exp(confint(modAccMaxPol)) # CI for paper
## Computing profile confidence intervals ...
##                   2.5 %      97.5 %
## .sig01         1.000000    1.310279
## .sig02         1.000000    1.646929
## .sigma         1.194421    1.392101
## (Intercept) 3308.135397 5250.567110
plot(modAccMaxPol)

# QQPlot for group-level effects
qqnorm(ranef(modAccMaxPol)$plant[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxPol)$plant[[1]])

# QQPlot for group-level effects
qqnorm(ranef(modAccMaxPol)$FlowerNumber[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxPol)$FlowerNumber[[1]]) 

# not reporting anthers, because they're basically the same as the pollen
modAccMaxAnth <- lmer(formula = log(a_T_anth) ~ (1|plant/FlowerNumber), data = md)
summary(modAccMaxAnth)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(a_T_anth) ~ (1 | plant/FlowerNumber)
##    Data: md
## 
## REML criterion at convergence: 19
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57600 -0.38059  0.05996  0.30044  2.00083 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  FlowerNumber:plant (Intercept) 0.01712  0.1308  
##  plant              (Intercept) 0.10281  0.3206  
##  Residual                       0.05554  0.2357  
## Number of obs: 32, groups:  FlowerNumber:plant, 15; plant, 7
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   8.2824     0.1355   61.12
exp(confint(modAccMaxAnth))
## Computing profile confidence intervals ...
##                   2.5 %      97.5 %
## .sig01         1.000000    1.408217
## .sig02         1.024878    1.843027
## .sigma         1.190678    1.405926
## (Intercept) 2988.212500 5284.547669
plot(modAccMaxAnth)

# diagnostics
plot(modAccMaxAnth)

# QQPlot for group-level effects
qqnorm(ranef(modAccMaxAnth)$plant[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxAnth)$plant[[1]])

# QQPlot for group-level effects
qqnorm(ranef(modAccMaxAnth)$FlowerNumber[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(modAccMaxAnth)$FlowerNumber[[1]]) #